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Method developed by Sandalov and coworkers [Int. J. Quant. Chem. 94, 113 (2003)] is applied 
to inelastic transport in the case of strong correlations on the molecule, which is relatively weakly 
coupled to contacts. Ability of the approach to deal with the transport in the language of many- 
body molecular states as well as take into account charge-specific normal modes and nonadiabatic 
couplings is stressed. We demonstrate capabilities of the technique within simple model calculations, 
and compare it to previously published approaches. 
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I. INTRODUCTION 

Development of experimental capabilities dealing with 
nanostructures brings necessity of appropriate theoret- 
ical description of quantum transport (charge, spin, 
and heat) in mesoscopic junctions to the forefront of 
research.— Indeed, a lot of work has been done in this 
direction. In particular, many approaches are based on 
the Landauer expression for current through such junc- 
tions in elastic tunneling regime.— One of specific features 
of molecular transport junctions, the focus of molecular 
electronics, is flexibility of the molecules, which results in 
inelastic features being much more pronounced in trans- 
port through such junctions as compared e.g. to semicon- 
ductor quantum dots. Inelastic features are used as a di- 
agnostic tool, helping to assure presence of the molecule 
and study its characteristics in the junction within in- 
elastic electron tunneling spectroscopy both in the off- 
resonance (IETS)2 and resonant tunneling (RIETS) 4 . sit- 
uations. Detailed discussion of the inelastic transport in 
molecular junctions can be found in Refs. 

Theoretical description of IETS is well established to- 
day both within simple models 7 -^ ' 10 ' 11 ' 12 and more real- 
istic calculations ! 13 ' 14 ' 15 ' 16 ' 17 ' 18 Ability to predict quan- 
titatively experimental findings is a sign of maturity 
of the field. From theoretical perspective this success 
is caused by ability to use well-established nonequilib- 
rium perturbation (in electron-vibration coupling) tech- 
nique. Indeed, in the off-resonant situation electron- 
vibration coupling M is an effectively small parameter, 
M < ^AE 2 + (r/2) 2 with AE being resonant off-set 
and r characterizing strength of molecule-contacts cou- 
pling. This allows expansion of the evolution operator in 
powers of M, and truncation at low (M 2 ) order, the Born 
approximation, is usually sufficient to get quantitatively 
correct predictions of IETS signal in molecular junctions. 

Resonant tunneling situation, AE = 0, provides richer 
physics. While weak electron-vibration coupling, M <C 
r, is treated within perturbation theory also here ) 19 ' 20 
the last fails in the opposite situation, M > T, when 
e.g. formation of polaron on the molecule becomes pos- 



sible. Theoretically this case till now mostly was treated 
either within scattering theory (or isolated molecule) 
approach ) 21 i 22 i 23 i 24 or within quasiclassical (rate or gen- 
eralized rate equations) scheme i 20 ' 25 ' 26 While the first 
treats electron-vibrational interaction (numerically) ex- 
actly, it disregards Fermi populations in the contacts, as 
well as dynamical features due to their presence^ 7 - and 
may lead to erroneous predictions ! 20 ' 28 The last disre- 
gards quantum correlations, and as such is applicable 
to either high temperature, fc^T 3> fku, (truly classical) 
situations, or (for generalized rate equations approach) 
quantum situations when correlations in the system die 
much quicker than electron transfer (between contact and 
molecule) time (~ 1/r). Besides, these schemes lack 
formal procedure for improvement of their results sim- 
ilar to taking into account higher order terms in per- 
turbative expansion. Recently we proposed nonequilib- 
rium equation-of-motion (EOM) approach perturbative 
in molecule-contact coupling capable of dealing with RI- 
ETS situation^ The approach is formulated for simple 
resonant level model, but can be easily generalized for 
more realistic situations. 29 While it incorporates contacts 
(and hence nonequilibrium character of the junction) into 
consideration, the price to pay is (generally) more ap- 
proximate level of description of electron-vibration in- 
teraction on the bridge. Alternatively, schemes explor- 
ing particular parameter regions (slow vibration loq < 
T^ZLZ 2 . small, V < wor 3 - or big, V > wor 4 - bias) were 
proposed at a model level. 

Another point especially important for resonant trans- 
port (both elastic and inelastic), when actual oxidation 
or reduction of the molecule takes place, is necessity to 
speak in the language of many-body molecular states con- 
trary to single-particle molecular orbitals (the last is used 
in most ab initio transport calculations today). This in- 
cludes electronic structure reorganization upon charging, 
state dependent vibrational modes, anharmonicities, and 
non-Born-Oppenheimer couplings. First schemes trying 
to treat transport in a many-body molecule states lan- 
guage were recently proposedj 25 ' 26 ' 35 

Difficulties in describing RIETS stem from absence of 



2 



well-established nonequilibrium atomic limit for molecule 
in the junction. Indeed, contacts play role of bound- 
ary conditions responsible for establishing nonequilib- 
rium state of the molecule. The last is a complicated 
mixture of different charge (and excitation) states. Ap- 
proaches developed in molecular electronics community 
so far either disregard boundary conditions and treat 
molecule as an equilibrium object (scattering theory and 
isolated molecule treatments), or establish this nonequi- 
librium state (mostly nonequilibrium Green function, 
NEGF, approaches) being unable to map it into separate 
charge (or more exactly state) constituents. Note, den- 
sity matrix based schemes capable of such mapping were 
developed recently ! 36 ' 37 Such schemes however, miss time 
correlations, which may become important in e.g. noise 
spectrum calculations. 

Hubbard operators is a natural language to talk about 
system (in our case subsystem - molecule) in terms of 
its states. Thus one seems to be interested in utilizing 
nonequilibrium Hubbard operator Green function tech- 
nique for description of situations similar to RIETS, 
when ability to establish nonequilibrium atomic limit 
of system is desirable. The approach should be ca- 
pable to provide a systematic way to take correlations 
into account (similar to perturbative expansion in stan- 
dard diagrammatic techniques). Such approach originat- 
ing from Kadanoff and Baym functional derivative EOM 
scheme, was developed in the form of equilibrium Hub- 
bard operator GFs by Sandalov and coworkers for ma- 
terials with strong electron correlations (magnets with 
localized and partly localized moments, Mott insula- 
tors, Kondo lattices, heavy fermion systems, and high-Tc 
superconductors)^ The method so far has been applied 
to model elastic transport through quantum dot a 39 ' 40 ' 41 
and lowest states of double quantum dots ; 42 ' 43 ' 44 ' 45 and 
is completely ignored in the molecular electronics com- 
munity. 

The goal of our present consideration is to introduce 
inelastic transport description in the Coulomb blockade 
regime within proper non-equilibrium atomic limit, and 
to attract attention of the molecular electronics com- 
munity to the proper nonequilibrium approach capa- 
ble of speaking in the language of many-particle states 
(rather than single-particle orbitals) and taking into ac- 
count both molecular charge state dependent normal 
modes (presently largely ignored in simulations) and non- 
adiabatic couplings. Note, that Kondo physics is be- 
yond the scope of current consideration. The approach 
takes into account only on-the-molecule correlations in a 
way that generalizes previous considerations ^ 5 ' 26 Note, 
that including many-body molecular states into consid- 
eration of transport potentially allows for: 1. Much more 
accurate molecular structure simulation out of equilib- 
rium than in current ab initio schemes, due to possibility 
to employ equilibrium quantum chemistry methods as 
a starting point for self-consistent procedure. 2. Proper 
treatment of oxidation/reduction and corresponding elec- 
tronic and vibrational molecular structure changes, as 



well as non-adiabatic couplings, 3. Ability to deal with 
general form of electron-vibration interaction as long as 
it is localized in space (in the spirit of Ref. but in 
addition retaining many-body character of the junction), 
4. Calculation of noise spectrum of the junction due to 
preserved time correlations (see e.g. Ref. |46| for detailed 
discussion), 5. Proper treatment of degenerate situations 
due to preserved space correlations (see e.g. Ref. for 
discussion). Structure of the paper is the following. In 
Section |TT] we briefly describe the method in terms of 
many-electronic states of the system, and compare it to 
previously proposed generalized master equation scheme. 
Section HTT1 presents numerical examples of its application 
to transport with discussion. Section HVl concludes. 



II. METHOD 

Here we introduce model of molecular junction, briefly 
review the basics of nonequilibrium Hubbard Green func- 
tion technique, and compare it to previously proposed 
generalized master equation approach. 



A. Model 

As usual we consider molecular junction consisting of 
3 parts: left (L) and right (R) contacts and the molecule 
(M). The contacts are assumed to be reservoirs of free 
electrons each at its own equilibrium. Molecule (or su- 
permolecule if inclusion of parts of contacts is required) 
is the nonequilibrium part of the system. Besides, any 
external potential, e.g. gate voltage probe, or additional 
contacts can be added to the picture if necessary. The 
Hamiltonian of the system is 

H = H l + H m + Hr + H t (1) 

where Hk (K = L,R) is Hamiltonian for contact K 

H K = ^2 £kc\'c k , (2) 

k£K 

Hm is Hamiltonian of the isolated molecule, and Ht is 
coupling between the subsystems 

H T = ^ \Vkmc\dm + VmkdlnCk) (3) 

k£{L,R};meM 

(v (d) and c 1 (c) are creation (annihilation) operators 
for electron on the molecule and in the contacts, respec- 
tively. Their indices m and k denote electronic state in 
some chosen single-particle basis, and incorporate all the 
necessary quantum indices (e.g. site and spin). 

Now we want to consider molecular subsystem in the 
basis of many-electron states \N,i >, where N stands for 
molecular charge (number of electrons or excess electrons 
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on the molecule) and i numerates different (e.g. excita- 
tion) states within the same charge state block. Gener- 
ally these states should not be orthonormal, and conse- 
quences of overlap between different molecular states (as 
well as overlap of molecular and contact states) were con- 
sidered in several papers . 38 ' 48 In what follows however we 
chose the states \N, i > to be orthonormal 



< N,i\N',i' >= S N N >5i. 



(4) 



in order to keep notation as simple as possible. Hubbard 
operators are introduced as usual 



X 



(N,i;N>,i') = \N,ixN',i'\ 



(5) 



In terms of these many-electron states transfer Hamil- 
tonian becomes 

Ht= Y (v k M~clx M +V Mk X^ M c k ) (6) 

ke{L,R};M 



where 



M = (N,i;N + l,j) 



(7) 



denotes transition of the system from state \N + l,j > 
to state \N, i >, while M = (N+l, j; N, i) stands for the 
backward transition. Transfer matrix element is 



V kM = Y Vkm < N >*\*™\ N + 1 >3> 



(8) 



and V Mk = V k * M . Often many-electron states are chosen 
as eigenstates of isolated molecule, in this case 

Hm = Y, En,iX( N ^ n ^ (9) 

\N,i> 

with Em,% energies of the isolated molecule states. 

B. Current expression 

Following derivation by Meir and Wingree n 49 ' 50 one 
gets usual expression for the current at interface K = 
L, R 

lR(t) =\ f d^Tr[ 

K J-oo 

S< (t, t x ) G> (ti ,t) + G> (t, h) E< (tj , t) (10) 
- £> (t, h) G< (h ,t)-G< (t, h) £> (h , t)] 

which for steady-state situation simplifies to 

"+°° dE , 



Ik = 



2n 



Tr[Z<(E)G>(E)-X>(E)G<(E)] 

(11) 

The only difference from the standard NEGF expres- 
sion is that Tr[. . .] in (flTJ)) and (TlTj) goes not over single- 
electron basis, but over basis of single- electron transitions 
A4, Eq.([7]), between many-particle states of the molecule. 



Self-energies Sk in (fTTJf and (jTTJ) are defined on the 
Keldysh contour as 

[^k(t,t')] mm , ee Y VMk9k{r,T')V k M' (12) 

keK 

with 

g k (T,T') = -i<T c c k (T)cl(r')> (13) 
GF for free electrons in the contacts. SEs projections are 

[Zk(E)] mm , = iT MM ,(E)f K (E) (14) 

[^k(E)] mm , = -i^MM'iE) [1 - f K (E)} (15) 
where 

Y MM ,(E) ee Y, v Mk VkM' S(E - s k ) (16) 

k£K 

and fn(E) is the Fermi distribution in contact K. 

GFs in dill) and ^TTJ) are Hubbard operator GFs de- 
fined on the Keldysh contour as 



G 



MM 



(r, t') ee -i < T c X M (r) X f M , (r') > (17) 



Note, that operators in (fT3"|) and (jTTJ) arc in Heisenberg 
representation. Note also, that M and M' in (fT7)) may 
be (in principle) arbitrarily far away from one another in 
the charge space. GF (jTTJ) represents correlation between 
different single-electron molecular many-body state tran- 
sitions due to coupling to the same bath (contacts). In 
practice however, it seems unreasonable to go beyond 
correlations between nearest charge space blocks. 

In order to show connection of the present formalism 
to previously proposed generalized rate equation (master 
equation in the Fock space) approach, we have to real- 
ize that the last misses correlations both in space and 
time. So to reduce present GF description to the mas- 
ter equation in the Fock space, we need to make several 
simplifications: 

1. Diagonal approximation. We have to stick to di- 
agonal elements of GF only Gmm with M. — 
(N,i;N+l,j). 

2. Markov approximation. We have to consider only 
GFs of equal times G(t, t). In order to reduce GFs 
of different times entering (jTTJ) to equal times quan- 
tities we use approximation 

Gmm (t - t') sa exp[i (t' - t)] G MM (t - t) (18) 



where 



- E 



N,i 



Now, noting that 



iG 



< 

MM 



t) = R 



(t-t) = R 



N 



N+l 



(19) 

(20) 
(21) 
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are probabilities to find molecule in state \N,i > and 
\N + l,j > respectively, we get from (fTTj) 



N;i,j 



(22) 



- r (AT,i;Ar+lj) [1 ~ fK(Ej - E i)]P. 



?N\] pJV+1 
j 



If now one restricts attention only to particular charge 
space block No and its nearest neighbors, one gets 
Eqs. (6) and (7) of Ref. [U 



C. General equation for GF 

Now, when expression for the current is established, 
we need a procedure to calculate Hubbard operators GF 
(fT7|) . Note, that standard diagrammatic techniques are 
inapplicable here, due to lack of the Wick's theorem 
(since Hubbard operators are many-particle operators). 
An alternative to diagrammatic expansion in the form of 
functional derivative equation-of-motion technique was 
developed in Ref. [H. Here we briefly review steps needed 
to obtain EOM for GF we will use in our numerical sim- 
ulations. 

Following Ref. [H we start by writing EOM for Hub- 
bard operator Xm(t), where M = (N,i;N+ 1, j). This 
leads to 



i— - A 



X m {t) 



f— Vk{N+l,j;N+2,l) c\{t) X ( N mn +2 ^){t) 

ke{L,R}-,e 

- Vk(N-i4-N,i) c fc (r)X(jv-i,^jv+i,j)(T) (23) 

+ V(N+lJ;N,e)k X(}f,i;N,t) i T ) Cfc(r) 
+V(N+l,t;N,i)k X(N+l,e-N+lj) (t) C k (t ) 



In what follows we disregard first 2 terms on the right- 
hand-side, since they describe simultaneous transfer of 2 
electrons between contact and molecule, which is beyond 
the scope of this consideration. It is clear that when 
writing EOM for GF (TIT)) terms in the right-hand-side of 
will produce correlation functions of the form 



<T c X 6 (r)c k (r)xl l ,(r')> 



(24) 



which can not be factorized into product of single- 
excitation GF (fTF)) and contact single-electron GF ([13")) 
due to lack of the Wick's theorem. 

In order to make this separation a trick with auxiliary 
fields (r) is employed. We need to introduce additional 
disturbance potential 

Hu(t) = ^ U(N,i;NJ)(t) X(Ar,. i; ]Vj)( T ) (25) 
N;i,j 

and corresponding generating functional 



Su = exp 



drHuir) 



(26) 



Then defining GF of 2 arbitrary operators A and B in 
the presence of auxiliary fields U 



G ab (t,t') =-i< T c A{t)B(t') > u 
= _.<T c S u A{r)B{r') > 



(27) 



<T C S U > 
one easily can get the following identity 

-i<T c X s (r")A(T)B(T')> u (28) 

S 



<T c Xs(t") > u +i- 



Gab(t, t') 



8U(:{t")_ 

Eq. ({2"8"|) allows to express correlation function in 
terms of single-excitation GF and its functional deriva- 
tives relative to auxiliary fields. Note, that putting (at 
the end) auxiliary fields to be zero turns (f2"7]) into a stan- 
dard definition of GF. 

So, introducing auxiliary fields as in Eqs. and (P2S1) 
and using expression (|28|) . one gets general EOM for Hub- 
bard operator GF (TiTl) in the form 



■ 9 . n 
l ~dr ~ 



Gmm'(t, t') 

~ E {M(N+l,j;N +!,{.) ( T ) G(N,i;N+l,t)M' ( r : T ') 



—M(N,£;N,i)( T ) G(jv^;AT+l,j) (t, t')) 

5(t,t')P mm ,(t) 



E 



< T c X(iv,i ; jv/) (t) >u +i- 



^ / dr" '£(N,t;N+l,j)M"{ T > T ") Gm"M'(t",t') 

s 



M 



< T c X( N+lt£ . N+ltj }(T) > u +i 



SU, 



(N+l,£;N+l,j 



X ^ / dT"Y l ( N i . N+1 i>) M n(T,T")GM"M'{T",T') J 
M" Jc ) 

where is defined in (I19| and 

PmM> =< T c X( Na . na ,){t) + X( N+1 .j'-N+l .j){ T ) >u 

(30) 

Eq. ((29|) is a general equation for Hubbard operator GF, 
representing an alternative to standard diagrammatic 
technique approaches. Role of expansion in small param- 
eter here play functional derivatives in auxiliary fields. 
Level of approximation is defined by order of the deriva- 
tive used in evaluation of GF. At the end of differenti- 
ations auxiliary fields are put to be zero, and resulting 
expression is equation for GF at particular level of ap- 
proximation. 



D. First loop approximation 

The simplest approximation, Hubbard I (HI), is ob- 
tained from (|29|) by keeping only diagonal averages, omit- 
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ting all functional derivatives, and hi — > 

m Gmm'{t,t') = 5(t : t')8mm> Pm (31) 

+ P M J2 l dT"Z MM „(T,T")G M „ M/ (T",T') 
M" Jc 

Following most of the papers, employed the method so 
far , 39 i 40 ' 41 in our consideration we go one step further - 
we take one functional derivative to get so called first loop 
approximation. Note, here we take derivative of the GF 
only, disregarding fluctuations of the spectral weight P. 
After performing the differentiation once more we keep 
only diagonal averages, omit all functional derivatives, 
and hi — > 0. Since the procedure was described in details 
in many papers (see e.g. Refs. HU12), here we present 
only final result 



i— - A 
Ot A ' 



M 



K ,e M" Jc 



[E(N,k;N+1J)M"{ t ) T ") D M"(N-l,t;N,i)(T" ,T +) 
x G(N-l,t;N,K)M'( T ) T ) 

- ^(N,K;N+t,j)M"( T i T ") P > M"{N,k-N+11){ t " \ T +) 

x G( N ^. n+1 ^ M i(t, t) 

+ ^(N.i:N+1,k)M"{ t i t ") Dm"{N,1;N+1,k){ t > r +) 
x G(Ar^ ; ]v+lj)X'( T > t) 

— ^{N,i;N+l,K,)M"( T i T ") D M" (N +l,j-N +2,1){ T " ■ T + 

x G(N+l, K ;N+2,t)M'{ T i T ')] 

= 5(t, t')8 M m> Pm 

^ / dr" EmM"(t,t")Gm"M'(t" ,t') 



(32) 



P 



M 

M" uc 

where M. = (N, i;N+l,j) and D is so called full locator, 
which (in the first loop approximation) obeys the same 
equation (|32|) as GF but without spectral weight Pm 
multiplying delta function in the right-hand-side. 

Expressions for GFs G and D (first loop approxima- 
tion) in the shorthand (matrix in both Fock space and 
Keldysh contour variables) notation can be written as 



b- 1 D 



P 
1 



where 



D~ 



Aa4 = 



i— - A M - PS 

OT 



A 



SA 



M 



(33) 
(34) 

(35) 
(36) 



and 8 Am is given by the second term in the left-hand- 
side of Eq. (|52")) . it is responsible for shifts of transition 
energies in the molecule due to contacts induced corre- 
lation. One sees that Eq. lfM]) has usual structure of the 
Dyson equation, which is obtained in standard diagram- 
matic expansion. The only difference is dressing of SE £ 



by spectral weight P. Thus formally one can use all the 
standard equations, using dressed SE S everywhere, to 
get desired projections of GF D. When D is known, G 
is obtained by simple matrix multiplication 



G = DP 



(37) 



Note, side of matrix dressing by spectral weight P is dif- 
ferent for E and G, compare Eqs. (|55|) and ([57)) . Note 
also, that the scheme is self-consistent, since both tran- 
sition energies shift 8 A and spectral weights P depend 
on GF G, while the last depends on these quantities. In 
particular (M = (N,i; N + 



Pm = N N j + N N+ ij 

NN,i =< -^(JV,i;JV,i) >= ^G MM (t,t) 



(38) 

for any j 
(39) 

Nn+ij ^< ^(N+i, r ,N+i,f) >= -iG MM (t, t) for any i 

(40) 

Here Ajv.i and Nn+ij are probabilities to find system in 
the state \N, i > and \N + > respectively. 



III. NUMERICAL RESULTS AND DISCUSSION 

Here we present results of simulations within first loop 
approximation. In order to speed up calculations we em- 
ployed also diagonal approximation^ We consider trans- 
port through quantum dot and double quantum dot, dis- 
cuss obtained data, and compare it to previously pub- 
lished results. 



A. Quantum dot 



In the case of quantum dot molecular Hamiltonian is 

(41) 



H M = e a h a + Uh^h l 

<7={U} 



where a indicates spin projection and n a — d\d a . Full 
Fock space of the molecular part of the system (without 
vibrations) consists of one empty state (|0 >= |0,0 >), 
two single-electron states, (| T>= |1,T> anQl I 1>= 
1,|>), and one doubly occupied state (|2 >= |2,0 >). 
Transitions between these states to be considered are spin 
up electron transfers (0 f and J. 2) and spin down electron 
transfers (0 J. and \ 2). Writing Eg. ([52"]) in the basis of 
these transitions one gets equations obtained in Ref. |40| . 

Figure [T] presents conductance map for elastic trans- 
port through quantum dot. Parameters of the calculation 
are T = 10 K, e a = -0.5 eV, = 0.01 eV (a =T, | and 
K = L,R), and U = 1 eV. As usually one has areas of 
blockaded transport (inside part of diamonds) with fixed 
population on the dot (0, 1, and 2 from right to left), and 
transition areas (between the diamonds) where popula- 
tion on the dot is noninteger (see e.g. Ref. [5l] for more 
detailed discussion) . 
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FIG. 1: (Color online) Conductance map for elastic transport 
through quantum dot. See text for parameters. 



In order to treat inelastic transport we add to Hm 
molecular vibration linearly coupled to electron(s) on the 
dot 



n, 



(42) 



Such model is frequently used to describe inelastic trans- 
port in molecular junctions. In a sense it is similar to 
Marcus theory, and describes shift of the molecular vi- 
bration when molecule is charged due to electron trans- 
fer from/to the contacts. In general, model with non- 
diagonal electron-vibration coupling also can be consid- 
ered within the formalism. 

After small polaron (Lang-Firsov or canonical) 
transformation 52 linear coupling term is eliminated, 
while energy level position e CT and Hubbard repulsion 
U are renormalized (e a — > e a — M 2 /ui n and U — » U — 
2M 2 /u>o), and transfer matrix elements in Ht, Eq.Q, 
are dressed with shift operators (d a — » d a X) 



X — exp [— A(a' — a)] 



(43) 



where A = M/luq. In what follows we disregard renor- 
malization of e a and U , assuming that it was included in 
definition of these parameters. 

Now molecule is characterized by direct product of 
electronic and vibrational spaces, so its state should 
be indicated by additional index v showing state of 
the vibration, i.e. molecular subspace is spanned by 
states \0,v >, | 1,v >, | l,v >, and \2,v >, where 
v E {0, 1, 2,3,.. .}. One has to consider the same elec- 
tronic transitions as in the case of elastic transport, 
but in addition all possible transitions between states 
of the vibration have to be included. Transitions be- 
tween these states (within the model) are possible only 
by electron transfer between molecule and contacts. Due 
to shift operators (|4"3"| appearing in Ht SEs (fT2")l are now 



dressed with corresponding vibrational overlap integrals 

(M = (N,i,Vi;N+l,j,Vj)) 



with 



< v\X\v' >= e- x2 / 2 (-l)( v - v 'W v - v 'h v " 



(44) 



(45) 



1/2 



where v m i n (v max ) is minimal (maximal) off and v' , 9{x) 
is step function, and L™ is Laguerre polynomial. Note 
an important formal difference between the present ap- 
proach and the one presented in Ref.[5l|. While in the last 
we had to consider separately electron and phonon dy- 
namics, which leads to convolution of electron GF (elec- 
tron dynamics) with Franck-Condon factors (phonon dy- 
namics), here the situation is different. Since we con- 
sider generalized Fock space (product of electronic and 
vibrational ones), within the formalism strictly speak- 
ing we do not have inelastic processes at all. Instead one 
has to consider elastic scattering events between electron- 
vibrational states. As a result the role played previously 
by the Franck-Condon factors (to introduce vibrational 
dynamics) now is included into Hubbard GF of the gen- 
eralized Fock space. 




FIG. 2: (Color online) Conductance map for inelastic trans- 
port through quantum dot. Linear coupling model (J42J) . See 
text for parameters. 

Figure [5] presents conductance map for inelastic trans- 
port through quantum dot within linear coupling model 
(|42l) . Parameters of the calculation are u>q — 0.2 eV 
and M = 0.4 eV, all the other parameters are as in 
Fig. [1] Within the calculation we restricted vibrational 
subspace to 4 lowest levels (v € {0,1,2,3}). As ex- 
pected, besides elastic peaks in the conductance map we 
get additional resonant vibrational features correspond- 
ing to inelastic processes. This Figure is equivalent to 
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Fig. 3a of Ref. |5j, where calculation was done within 
perturbative (in coupling to electrodes) nonequilibrium 
EOM approach for the same model. As previously,— 
within the model (see also discussion below) distance 
between the diamond edges (elastic peak) and vibra- 
tional sidebands is defined by the oscillator frequency. 
Increase in electron- vibrational coupling would result in 
both more pronounced vibrational features and suppres- 
sion of transport in the low source-drain voltage region 
due to Franck-Condon blockade. While increase in tem- 
perature would produce also vibrational sidebands corre- 
sponding to phonon absorption (features inside the dia- 
mond). 

Finally, we want to demonstrate capabilities of the 
present scheme, which go beyond approaches previously 
used to treat inelastic transport. Suppose our molecule is 
small enough, so that upon charging it changes its normal 
modes essentially. Suppose also, that from all the normal 
modes of the molecule only one is coupled to tunneling 
electron. Inelastic transport in this case can be modeled 
by assigning different vibration frequencies to different 
charge states of the molecule. In our quantum dot model 
this corresponds to situation, when vibrational frequen- 
cies for |0, v >, \a,v >, and |2, v > states are different - 

uJq , u>q X \ and u>^ respectively. Self-energies due to elec- 
tron transfer between molecule and contacts once more 
have to be dressed by overlap integrals between different 
vibrational wavefunctions, as is shown in Eq. fH]) . How- 
ever this time (when vibrational frequencies change) the 
integrals should be calculated in the way discussed in 
Refs. [HH 




FIG. 3: (Color online) Conductance map for inelastic trans- 
port through quantum dot. Charge state dependent frequen- 
cies. See text for parameters. 

Figure [3] presents conductance map for inelastic trans- 
port through quantum dot, when vibrational frequency 
depends on charge state of the dot. Parameters of the cal- 



culation are u>^ = 0.2 eV, uj^ 
and shift vector for both transitions was taken to be 



0.3 eV, 



0.25 eV 



0.5 A(see Refs. I531I541 for detailed explanation), all the 
other parameters are as in Fig. [T] One sees that result of 
calculation is counter intuitive at the first sight. Naively 
one could expect to see inelastic peaks at each diamond 
edge (each charge state of the quantum dot) being sep- 
arated by frequency corresponding to the neighboring 
charge state (RIETS probes frequencies of the interme- 
diate ion). Real picture is more complicated however. 
Let consider electron transfer between 2 particular charge 
states of the quantum dot, say between states |0, vq > and 
|er, v\ > upon electron transfer from contact to molecule. 
In this case change in the subsystem energy, which will be 
observed in transport as inelastic peak in conductance, is 

v\uj^ — vquj^ (we omit here change in elastic electronic 
energy for simplicity, this will define only position of the 
elastic peak in conductance). Since Vq and v\ in princi- 
ple can be any non-negative numbers, it is clear that one 
can observe a a progression of frequencies. Note, that 
in this progression one can see inelastic peaks in conduc- 
tance, separated from the elastic one by frequency which 



does not exist in the system at all (e.g. 



')• 



Note also, that due to overlap factors involved the lowest 
frequencies of the progression will be observed better in 
RIETS signal. Non-adiabatic couplings can be included 
in calculation in a similar way. 



B. Double quantum dot 

Molecular Hamiltonian for double quantum dot is 
H M = ^2 e la h la - t 12i<T [d\ a d 2a + 4<r^ 

i={l,2}ff={T,l} 

(46) 

+ 2J U l h^h i i + Ui 2 nifi2 
i 

where i — {1,2} numbers sites and a = {f , 1} stands for 
spin projection, <v (d) is creation (annihilation) operator, 



d]„dia, and n. 



HI 



Hi 



We assume that site 1 



is coupled to the left contact, while site 2 - to the right. 

We chose many-body states for molecular subsystem 
in the form ll j, 1 J.,2 |, 2 |>. Unlike choice of 

Refs. smug! these are not eigenstates of the molecu- 
lar Hamiltonian. As a result EOM for Hubbard operator 
GFs couples them also by hopping t\2,v Besides this all 
the treatment presented in Section HT1 remains the same. 
There are 16 states (1, 4, 6, 4, and 1 states for 0, 1, 
2, 3, and 4 electrons in the system respectively) and 32 
single-electron transitions (16 for each spin block) to be 
considered. 

Figure 2] shows conductance map for clastic transport 
through double quantum dot. Parameters of the calcu- 
lation are T = 10 K, e il7 = -0.5 eV, t 12 a = 0.01 eV, 

rf ff = r* = o.oi ev, r* = r& = o,u^=u 2 = u = 

1 eV, U\2 = 0.5 eV, and Ep = 0.5 eV. As usually one 
sees pattern of blockaded and allowed transport regions. 
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IV. CONCLUSION 




-0.5 0.5 

v„/u 



FIG. 4: (Color online) Conductance map for elastic transport 
through double quantum dot. See text for parameters. 



However, here this pattern is more complicated than in 
the case of quantum dot. Figure[H]demonstrates this pat- 
tern for a horizontal cut of Fig. [4] at source-drain voltage 
Vsd/U = 0.25. Shown are current (a) and probabilities 
to find the molecular subsystem in different occupation 
states. 




FIG. 5: (Color online) Elastic transport quantum dot at fixed 
source-drain voltage V s d/U = 0.25: (a) Current vs. gate volt- 
age, (b) Probability to find double quantum dot empty (dot- 
ted line, black), singly- (dashed line, green), doubly- (solid 
line, red), triply- (dash-dotted line, blue), or fully-occupied 
(dash-double-dotted line, magenta) . Parameters are the same 
as in Fig. 0] 



Approaches based on renormalization group technique 
(NRG, DMRG, real time RG, etc.) are routinely used to 
treat quantum impurity systems, where e.g. correlations 
between localized impurity (quantum dot or molecule) 
and delocalized states (contacts) leads to Kondo effect. 
Mostly these approaches were applied to description of 
strongly correlated systems in equilibrium. Descrip- 
tion of transport is more problematic, since one needs 
to know spectral function at finite temperature, where 
all excitations may contributed Nevertheless first ap- 
proaches dealing with transport started to appear as 
well j36 |55,i56 |51j58 i59 The main complication with imple- 
mentation of these methods (besides DMRG) for ab initio 
calculation is their complexity, so that all the calcula- 
tions done so far are restricted to simple models only. 
In the case of molecular junctions most of ab initio cal- 
culations done today are performed at the mean mean 
field level of treatment with effective single particle or- 
bitals used in place of molecular states. Such approach 
clearly breaks down in the resonance tunneling regime, 
where actual reduction/oxidation of the molecule lead- 
ing to corresponding electronic and vibrational struc- 
ture change become possible. Necessity of treating this 
regime in the language of many-body molecular states, 
thus incorporating on-the-molecule correlations, was re- 
alized and first approaches like e.g. generalized master 
equation approac h 25 i 26 where proposed. Here we gen- 
eralize this consideration by incorporating many-body 
molecular states language into nonequilibrium Green's 
function framework. The main formal problem here is 
that many-body states language makes the Wick's theo- 
rem inapplicable, and thus standard nonequilibrium di- 
agrammatic techniques can not be used. A workaround 
based on functional derivative equation-of-motion tech- 
nique for Hubbard operator GFs was developed by San- 
dalov and coworkers 3 ^ for equilibrium case. The method 
so far has been applied to model elastic transport through 
quantum dot a 39 ' 40 i 41 and lowest states of double quantum 
dots ; 42 ' 43 ! 44 ' 45 and is completely ignored in the molecu- 
lar electronics community. Here we employ the approach 
to deal with inelastic transport through molecular junc- 
tions in nonequilibrium atomic limit. We formulate the 
method within basis of charged states of the molecule. 
We demonstrate its ability to deal with transport situ- 
ation in the language of these states (rather than effec- 
tive single-electron orbitals), as well as take into account 
charge-specific normal modes as well as nonadiabatic cou- 
plings. Capabilities of the technique are illustrated with 
simple model calculations of transport through quantum 
dot and double quantum dot. Extension to realistic cal- 
culations is the goal of our future research. 
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Note, that approximation like the one presented in Sec- 
tioning] may result in unphysical behavior. In particular, 
retarded and advanced SEs and GFs are not Hermitian 
conjugates of one another. While the issue does not arise in 
the diagonal approximation implemented for calculations 
here, this is a problem for a more general consideration. 
A simple workaround is to use average of two Dyson-like 
expressions: one with -D -1 , Eq. (|35[) , applied from the left 
and one from the right. This leads to set of equations for 
GFs, which under Markov approximation reduce to widely 
employed equations for DM in dissipative environment. 



